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Abstract 

Sphere packings are essential to the development of physical models for powders, com- 
posite materials, and the atomic structure of the liquid state. There is a strong scientific 
need to be able to assess the fit of packing models to data, but this is complicated by the lack 
of formal probabilistic models for packings. Without formal models, simulation algorithms 
and collections of physical objects must be used as models. Identification of common as- 
pects of different realizations of the same packing process requires the use of new descriptive 
statistics, many of which have yet to be developed. Model assessment will require the use of 
large samples of independent and identically distributed realizations, rather than the large 
single stationary realizations found in conventional spatial statistics. The development of 
procedures for model assessment will resemble the development of thermodynamic models, 
and will be based on much exploration and experimentation rather than on extensions of 
established statistical methods. 

1 Introduction 

A disordered sphere packing process can be loosely defined as the type of spatial stochastic 
process whose realizations describe the final positions of a collection of spherical objects 
which are tossed into a container and come to rest. Sphere packings are used in physics 
and engineering to model the atomic structure of amorphous solids, the internal structure of 
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composite materials, and the static and dynamic properties of powders. They are examples 
of spatial stochastic processes which need to be fit to data in order to be properly used, yet 
no formal procedure has yet been developed for objectively assessing the quality of fitted 
models. These processes have a very different within-realization disorder to that found in 
processes normally studied in spatial statistics, and so very little study of these processes has 
been undertaken by statisticians or probabilists. Most work on packings has been undertaken 
by physicists and engineers who tend to focus on the mean behaviour of the processes and to 
exclude consideration of their variability, although their work has highlighted a serious need 
for fit assessment. By combining ideas and methods from spatial statistics and from physics, 
it is possible to develop a formal approach to inference which is useful in a scientific and 
engineering context. This approach to inference is also useful in the context of developing 
methods for fitting spatial models to highly dependent spatial data which is not amenable 
to fitting by conventional statistical methods. 

It will be generally be assumed that the spheres in packings are of equal diameter 
(monodisperse) rather than being of many different sizes (polydisperse). Spheres in M 2 will 
be termed discs, spheres in M 3 will be termed spheres, and in higher dimensions they will 
be termed hyperspheres. In the applied literature, planar packings of discs are sometimes 
referred to as packings of rods are cylinders. This is misleading, since it assumes absolute 
rigidity for the physical rods and cylinders being modeled. A physical packing is a collec- 
tion of spherical objects assembled by tossing those objects into a container and allowing 
a force to form them into a jammed static arrangement. Packings will be assumed to be 
disordered and not based on a point lattice unless otherwise stated. Ordered packings with 
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point lattice structure are of use in the study of geometry in high-dimensional spaces (111; |2j; |3|) 
and in the development of efficient coding (4), but these applications have no relevance to 
the modeling of disordered materials. Packings are also of great interest in more abstract 
mathematical contexts (e.g. the theory of discrete analytical functions (5)), but these uses 
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have no applications to materials either. 



1.1 Examples of use 



Sphere packings were first used in attempts to explain thepatterns seen in X-ray diffraction 
studies of liquids and glasses ((J 7). Bernal (8) and Scott (9) attempted to model macroscop- 
ically what an amorphous mass of densely arranged atoms might look l ike, using physical 
packings of monosized steel spheres as their model. More recent work (jlOl ) has suggested 
that this approach to modeling may not be feasible, since the packings of spheres are formed 
by the action of a different set of forces from those found at the atomic scale. If physical 
models can be shown to explain the patterns seen, then this use of the physical model would 
be analogous to the use of linear models in other sciences where the use of a simple model 
can be useful, even if it is partly incorrect. 

Sphere packings are more commonly used in the modeling of granular materials at rest. 
Granular materials are collections of solid particles surrounded by fluid or gas which form 
packings at rest, but which flow like liquids or gases when a sufficiently large force is applied. 



When put into motion, the complex interaction between the grains makes it impossib 
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model granular flow by means of conventional methods from statistical physics (I 111 ; 
Packings are essential components of models for the static and dynamic behaviour of granular 
materials which are used the study of avalanches, of mudslides, of soil liquification and 
erosion, and of the flow of powders through pipes and hoppers. 

Sphere packings can be used in the study of composite materials. In the study of metal 



sintering (1l4l). a packing can represent the initial state of a metal powder before it is com- 
pressed and heated to form a solid metal part. Polydisperse packings can be used as models 
for some some colloids and for concrete. A colloid is a mixture of two immiscible liquids, in 
which one liquid may tend to form spherical inclusions within the other. For some compo- 

rn 

sitions, the spherical inclusions in a colloid may appear to form a packing (1l5l ). Concrete 



is a composite of rocks of many different sizes held together by a matrix of cement (]16j). 



yidealized, but may be the only feasible 



Sphere-packing-based models for concrete are high 
way to abstractly model such a complex material (1171 ). 

Packings can be used as models for porous structures. A packed bed reactor is a large 
steel vessel filled with solid catalyst-impregnated particles. Liquid reactants flow in at one 
end of the vessel, and the product emerges from the opposite end. Design of these reactors 
requires being able to model both the reaction in the pore spaces around the particles and 
the transfer of heat through the particles. Sphere packings are the simplest form of particle 



packing that can be used to model the internal structure of these reactors (Il8l ) 



2 Defining a Packing Process 

A sphere packing process is a stochastic process whose realizations are arrangements of 
spheres which closely resemble the arrangements of spherical objects or atoms produced by 
physical processes in nature. Ideally, it would be possible to define the process mathemat- 
ically and then to analyze packing processes through mathematical arguments, as can be 
done for Poisson processes and spatial Markov processes. Insurmountable difficulties posed 
in defining the packed state and in defining the process make it impossible to undertake an 
analytical approach to inference. 

The first problem is to define what it means for spheres in a realization to be packed. 
Consider a configuration of spheres in M. d . Construct the Delaunay triangulation of the sphere 



centres (Il9l ). An interior sphere is any sphere whose centre has no triangulation edges linked 
to any of the synthetic points used to construct the triangulation. Boundary spheres are 
spheres which are not interior spheres. Identify all points of contact between spheres and their 
neighbours, and define the contact network to be the subset of the triangulation which joins 
centres of spheres in contact. If confining surfaces are present, then the contact network must 
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be augmented by edges connecting points of confining surface contact with the corresponding 
sphere centres. A sphere is jammed if it cannot be translated without overlapping another 
sphere. A boundary sphere is jammed if it cannot be translated without overlapping another 
sphere or penetrating a confining surface. Otherwise, it is a free boundary sphere. 

Given an configuration of spheres, then the simplest definition of a packing would require 
that all interior spheres be jammed and that the contact network consist of a single connected 
component. This definition is too strict, in that it is possible to have boundary spheres which 
are not jammed (called rattlers) inside of cages of jammed spheres which may also include 
confining surfaces. The simple definition is also not strict enough in many applications, in 
that it may include arrangements of spheres which are unstable when the packing is under 
the influence of forces. If forces are acting on the packing, then all of the spheres including 
those on the free boundary need to be in mechanically stable positions. 

Given these definitions, packings of monosized spheres can only be shown to exist by 
mathematical arguments if the spheres are arranged in a regular lattice structure, such 
as a hexagonal or cubic close packing. There is no way to demonstrate mathematically 
the existence of a disordered packing of monosized spheres. Spherical physical objects can 
be packed by mechanical processes, but the physical objects involved are neither exactly 
monosized or exactly spherical. Given any list of sphere centre locations in a possible packing, 
the exact distance between the points can never be exactly known, and so it cannot be said 
for certain if the list represents a disordered packing of monosized spheres. It is possible that 
there is no such thing as a disordered packing of monosized spheres, which would complicate 
the formal definition of a packing process. The same difficulties with existence apply to 
polydisperse packings, with some exceptions related to circle packings in the plane p). 

If disordered packings could be shown to exist, then it would be necessary to formally 
define a stochastic process which had disordered packings as its realizations. One candidate 



for such a model is to use the Gibbs process (1201 ) to model the sphere centre positions, which 



5 



would require that the packing be of infinite extent and be stationary. Gibbs processes are 
only useful when the only interactions between spheres are pairwise (12 ll ). In a Gibbs model 
for a packing, there would be many higher-order interactions which could neither be ignored 
nor simply approximated. 

Instead of formally defining a Gibbs model for a packing process, it may be possible 
to approximate it by explicitly listing the possible configurations for the spheres and then 
assigning a probability that each configuration will appear in a realization of a packing 
process. The construction of such configurations would require arguments similar to those 
used in the identification of the largest fraction of a square (or equilateral triangle or regular 
hexagon) which can be covered by non-overlapping monosized discs. These calculations 
can be done analytically for 9 or fewer discs, and computer-based proofs are needed for 10 
I). For larger configurations, approximate simulation-based arguments are 



to 27 spheres ( 



required ((23 



22 



2J) which makes this approach futile for any reasonably-sized packing. 
The lack of a useful formally defined model such as the Poisson process or the Gibbs 
process implies that likelihood-based inference is impossible for any physical packing process. 
Any inference or model assessment procedure must be based on replications of physical 
processes or on independent realizations generated by computer programs. 



3 Useful Models for Packings 

Since mechanically stable packings of physical objects exist in nature, it is scientifically 
sensible to attempt to represent them by some type of abstract model which is simpler 
to deal with than the physical objects themselves. While formal mathematical modeling is 
impossible, other physical systems and complex computer programs can be used as stochastic 
models instead. 
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3.1 Physical Models for Packings 

A physical model for a packing of objects is a different physical system which is easier to 
replicate and study than the original packing. Examples include using packings of steel 
spheres to represent the positions of atoms in amorphous solids or to represent the internal 
structure of catalyst pellets in a packed bed reactor. In the latter case, a scale model of the 
real system could be prepared out of different materials from the reactor, so that internal 
images of packed bed could be found. Physical models are also useful for investigating 



packings of soft objects and objects of irregular shape (1251 ) when these objects are too 
difficult to simulate. 

The physical model is stochastic because its realizations are generated by a process that 
can be modeled by a dynamical system which is sensitive to initial and boundary conditions. 
Pouring spheres from one container to another is an example of this type of physical process, 
since repetitions of the process under the same general experimental conditions will produce 
different disordered packed configurations. For any such process, there is a sample space of 
possible packed object configurations and a density which represents the relative likelihood 
of each configuration. The goal of model selection in the case of the physical model is 
to match the sample space and distribution to that of the physical system which is to be 
modeled. Both of these sample spaces and densities are inaccessible to theoretical argument 
and difficult to approximate, and so any inference must be based on comparisons of the 
replications. 

As involved attempts to model the molecular structure 
8|) both used packings of monosized steel spheres as 



Tlie first major use of physical mode 



of liquids (6>). Scott (J9) and Bernal (|7|; 
models for atomic structure. Packings were made inside of thick rubber balloons bound with 
rubber bands in order to maximize packing density and reduce the effects of gravity. Bernal 
also filled the balloons with black paint, which was then allowed to set. The agglomerated 
mass was dissected carefully in order to identify which spheres were in contact, and physical 



models of the contact network were constructed. More elaborate measurement methods were 



developed (1261 ). but interest in this approach soon waned on account of the amount of time 
and effort required to obtain useful results. 

Physical packings are also distinguished by the material that fills the voids around the 
spheres. A physical process that packs spheres in room temperature air is distinct from one 
that packs spheres surrounded by warm liquid wax or paint. To make low-density packings, 



spheres have been packed in a fluid of equivalent density (ll5l) . It also possible to make 



something very much like a packing by combining two immiscible fluids, one of which forms 



near-spherical inclusions in the other (1271 ) 



Physical packing processes for discs can also be developed. These could be as simple as 
swirling a tea tray full of thick coins and then inclining the tray until the coins come to 
rest. Processes of this kind produce realizations which are very close to hexagonal packings, 
and so more complicated methods are needed to obtain more disordered realizations. Non- 
overlapping discs can be placed at random locations on a slightly stretched sheet of rubber. 
The tension in the sheet is released, forcing the discs together into a smaller area. The 
positions of the discs are recorded, the sheet is retensioned, the discs are placed back on the 
sheet, and the tension is released again. By cycling through this procedure, a disc packing 



can be generated (1281 ). 



Interest in physical packings has been revived as new remote sensing methods have been 



dev elop ed. X-ray tomography can be used to map the interiors of very large packings ((29 



30 



31 



32j). Magnetic resonance imaging (1331 ) may be less expensive, but it also produces 
lower-resolution images and cannot be used on very large packings. Confocal microscopy has 



been used to study both sphere packings ([34 



35l ) and colloidal packings ( 1271 ) 
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3.1.1 What can be learned from physical packings 

Any model of a physical process must obey the constraints on packed configurations which 
have been revealed by the study of packings of physical objects. 

The choice of experimental procedure which defines a physical packing process has a 
great impact on the sample space of configurations associated with that procedure. Experi- 
ments with monosized spheres poured into containers showed that the density of the initial 
configuration of the pcn r cou.d be made ,nore dense by redding the packing or by vibrating 
it up and down (J9_|). There appeared to be a well-defined limit to the density achievable by 
these methods, which resulted in the proposition that there was a well-defined physical state 



termed a dense random packing of spheres. This state has been shown to not 



3e well-defined 



(1361 ). and it is possible to use mixing procedures in a Couette- Taylor cell (1371 ) and other 



methods (I38l ) to produce denser near-ordered structure packed configurations. The implica- 
tion for modeling is that an arbitrary packing model cannot be assumed to be sampling from 
the same sample space of configurations as does the physical process that is being modeled. 

The apparent existence of a well-defined dense random packing opened the question as 
to whether or not a disordered packing could be denser than a hexagonal close packing 
of spheres. The mean volume fraction of dense random sphere packings was thought to 
be between 0.63 and 0.64, but that did not rule out the existence of some configurations 
which might have higher volume fractions. For dimensions 2 < d < 8, the densest lattice 
packings can be identified. For discs in the plane, the densest packing of any kind can be 
analytically proven to be the hexagonal close packing (39). A computer-based proof with the 
same conclusion has been constructed for spheres (J40j). In dimensions 10 and higher, some 
non- lattice packings have been found which are denser than any known lattice packings (jj). 

Polydisperse packings are far more disordered than monodisperse packings. One way of 
disrupting the near-order of physical packings of monosized discs is to randomly introduce 
a small amount of size variation into the discs. Physical packings of a fixed sphere size 
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distribution can be very difficult to assemble, since the smaller spheres often segregate and 
accumulate at the base of the specimen. These difficulties also occur in simulated polydis- 
perse sphere packings (14 ll ; |42| ) . 



3.2 General Comments on Simulation Models 

A simulation model is a program which takes output from a random number generator and 
produces an arrangement of spheres which is close to that of a packing. While these algo- 
rithms could in theory be represented by a formal probabilistic model, they are sufficiently 
complex that no such model could be stated or usefully analyzed. 

3.2.1 Machine Precision and the Contact Network 

Any simulated monodisperse sphere packing will be defined by a list of sphere centre loca- 
tions. With the exception of some cubic packings, most of these numbers will be irrational. 
In any simulation, sphere location coordinates can only be specified to a finite number of 
decimal places. It is not possible to determine if the list of points produced by any simulation 
program form a packing or not. It is generally assumed that any simulated packed configu- 
ration can be collapsed to a rigid jammed configuration, but it is not possible to determine 
exactly which spheres in the collapsed configuration are in contact and which ones are not. 
Examination of physical packings shows that they contain many pairs of spheres which are 



close but not in contact (la: l43l) . If the simulation is intended to study a physical phenomenon 
which is very sensitive to the structure of the contact network, then this uncertainly could 
contribute additional uncertainty to any physical properties estimated from the simulated 
packing. 

If round-off errors are the only type of error present in the simulation, they can be min- 
imized through use of maximum precision variables in the code. It is also possible to adapt 
methods used to deal with errors arising from detectors and from Fourier reconstruction in 

10 



X-ray tomographic imaging of physical packings. In those reconstructions, the estimated 



distance between the closest neighbours was found to be approximately Gaussian (129l ). This 
distribution was used to estimate the mean number of neighbours in contact, and could be 
used to define a stochastic or deterministic method for deciding if two spheres are in contact. 

3.2.2 Boundaries 

Most physical packings are formed by a confining surface. This may be a plane orthogonal 
to the direction of the force driving the physical flow that produced the packing, or it may 
be a completely closed surface that jams every sphere within the packing. Simulations often 
omit confining surfaces, either to simulate small parts of much larger packings or to avoid 
difficulties in coding. 

It is possible to avoid using any boundaries. Spheres can be packed on the surface of 
a hypersphere of one dimension higher (144 ). but this can introduce unwanted effects from 
hypersphere surface curvature unless many spheres are used. Packings can also be assembled 
by simulating a force acting towards a central point. It is not clear that a packing simulated 
by a central force would represent any packing found in nature. 

Confining surfaces can be avoided by imposing periodic boundary conditions. Consider 
a packing in the plane and a rectangle W. If the discs are monosized with diameter 5, 
then any disc whose centre lies within 5/2 of a side of W has the part of the disc outside 
of W appear on the opposite side of the rectangle. Since this can occur on any edge, the 
window W is toroidal and the plane is tiled by copies of W containing a disc packing. If the 
boundary effects vanish after some distance k8 from the periodic boundary, then the interior 
of a realization may be indistinguishable from that of a stationary and isotropic packing. 



49 ). but these are often based 



Efforts have been made to determine a useful value for k (145 : 
on a single statistic and may not reflect every effect of periodic boundaries on the structure 
of a realization. 
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Periodic boundary conditions may be used within confining boundaries in order to save 



computational 



effort. Within cylindrical confining boundaries, it is possible to pack four 



periodic cells (1471). In any small system, this symmetrical structure may have a strong effect 



on the distribution of statistics calculated from the packing. 
3.2.3 Packings of non-spherical objects 

Most simulation models pack monodisperse spheres, since the code for these programs is 
relatively easy to write. It is dangerous to assume that such models can represent packings 
of more complicated objects, unless statistical methods can be used to show that simpler 
models can represent the physical phenomenon of interest. 

Modeling packings of randomly shaped objects is very difficult. While a sphere requires 
only its centre and radius to completely specify it, a reasonable simulation of a small rock 
might require a long list of numbers to specify the polynomial surfaces that are spliced 
together to produce a simplified description. When monodisperse spheres are packed, feasible 
arrangements can be determined by comparing between-centre distances to the diameter of 
the spheres. For random shapes, very complex code is needed in order to efficiently determine 
if any overlaps have occurred. Until these problems can be overcome, the best results for 
random shapes will come from studying physical packings rather than simulated ones. 

It is possible to simulate the packings of simple non-spherical objects. Ellipses can be 
packed in the plane and ellipsoids can packed in space because these shapes can be easily 
described and overlaps can be easily determined. Ellipsoids can be shown to pack more 



densely than spheres in simulations ( 148 ; 



49 



50l ) and in physical experiments (15 ll ). The 



shapes of ceramic rings have been approximated by assemblies of triangles ( 1521 ) as part of a 
simulation of packed beds. 
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3.3 Ballistic Algorithms 



Ballistic algorithms are the simplest algorithms which can be used to produce configurations 
of packed spheres or discs. They were the first developed and are relatively simple to code. 
The algorithms have no basis in physics, and packings generated by these algorithms are 



less dense than those generated by physical packing processes (1531). Ballistic models must 
be used with caution in physical modeling, as there is no reason to believe that they have 
the same sample spaces as do physical processes. 



55l ) to model the formation 



The first ballistic algorithms were developed by Void (154j ; 
of floes. Spheres are dropped at random locations onto a surface. They stop falling either 
when they hit that surface or another sphere. If they hit another sphere, then they either are 
locked in place with probability p, or are rolled down the packed spheres until they reach a 
stable position. With p = 0, the algorithm can produce realizations of loose packings. Void 



also developed a central version of this algorithm (1561 ) and one involving rods composed of 



a set of k spheres in contact with their centres arranged along a line (15 71 ) . 



The first high-density central packing simulators built up packings by placing spheres at 



59T ). These algorithms produce 



the closest possible position to the centre of a packing (158 ; 
packings with estimated volume fractions between 0.61 and 0.62. For central disc-packing 
algorithms, near-ordered patterns can be avoided by seeding the realization with a non- 



triangular configuration of discs (j60l ). 



Visscher and Bolsterli (16 ll ) developed a non-central algorithm which added periodic 
boundary conditions perpendicular to the base. A sphere is dropped at a random loca- 
tion, and falls vertically until it contacts another sphere. Then, it follows the shortest path 
along the surface of the already-placed spheres until it comes to rest in a gravitationally 
stable position. Each sphere drop is repeated k times, and the final position chosen is the 
position that is closest to the base. The algorithm produces realizations with estimated 
volume fraction of 0.582 in M 3 , which is less dense than a settled physical packing. In the 
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planar version of this algorithm, the initial spheres dropped have a slightly different size in 
order to prevent the formation of near-ordered realizations. 

The Visscher-Bolsterli algorithm has been extended in various ways. A central version 



has been developed (1621 ). with the intention of producing packings that contain bridging 



structures. Non-periodic boundaries have been added to the algorithm, so that the spheres 



fill a box in space (1631 ) or in the plane (1641), or fil 



boundaries, irregular shapes can be packed (J52 



68|). 



a cylinder ((65 



66 



671 ). Using periodic 



The Visscher-Bolsterli algorithm can be modified to pack spheres of different sizes. When 
the largest and smallest spheres are greatly different in size, the smaller spheres tend to 
accumulate at the base of the packing unless they are allowed to randomly stick at unstable 



positions (14 ll ). More complicated deposition rules have been devised to avoid size segregation 



in the plane ((69 



3) 



The density of physical packings can be increased by shaking them. The density of 
ballistic packings can be increased by subjecting the initial packing to rearrangements which 
superficially imitate shaking. A simple rearrangement algorithm orders the initially packed 
spheres by height, then r edep osits them at the same planar location using the Visscher- 



Bolsterli placement rules ((71 



721 ) . This algorithm was used to study segregation by size 
after shaking. In a shaking algorithm for monosized spheres, each sphere in the packing is 
displaced upwards by a small Normally distributed perturbation, and is then subjected to 
many small three-dimensional Normal perturbations which are allowed if no collisions occur. 
Once the number of collisions reaches a set threshold, the packing is collapsed from the 
bottom using VB deposition rules. This algorithm can increase the estimated mean volume 
fraction from 0.581 to 0.590, and still results in a loose packing ^7^ ). If the initial vertical 
perturbations are not sufficiently small, then the rearrangement can simplify the contact 
network and decrease the volume fraction. Rearrangement methods can also be applied to 



packings of irregularly shaped objects in a container with hard boundaries (174J). 



14 



3.4 Rearrangement Algorithms 

Algorithms which rearrange the points in point patterns were the first to achieve volume 
fractions similar to those seen in physical packings. These programs begin with a realization 
of a Poisson or a regular point pattern. The points in the pattern are subjected to determin- 
istic or random translations which eventually produce a sphere configuration close to that 
of a packing. These methods do not attempt to replicate the dynamic interactions between 
spheres which occur during the formation of a physical packing. 



The Jodrey-Tory algorithm ((75 



76 



771 ) was the first simulation algorithm to produce 
realizations with estimated volume fractions similar to those for dense physical packings. 
It is initialized by a configuration of n points in the interior of a rectangular prism with 
periodic boundaries. A sphere of unit radius is attached to each point. In the first stage of 
the process, the radius of each sphere shrinks by 0.0001 units. The distance between each 
sphere and its nearest neighbour is found, and the closest pair defining overlapping spheres 
is identified. This pair are moved apart along a line joining the points until the spheres no 
longer overlap. When the distance between the closest two spheres drops below a threshold, 
all spheres have their diameter increased by 0.0002 units, and the process repeats. After 2000 
cycles, a second routine of shrinking and translation removes all remaining overlaps. Once 
the initial configuration of points has been chosen, the algorithm is entirely deterministic. 



There are also versions of the algorithm which pack ellipsoids (781 ) and discs of two different 



sizes (I79l ). 

The Jodrey-Tory algorithm can be initialized with one of its own previous outcomes. If 
this done many times, then the configuration crystallizes to become denser and less disordered 



(l80l). This phenomena is not seen in the simplest physical processes for generating dense 
physical packings, and it suggests that the Jodrey-Tory algorithm is not sampling from the 
set all feasible packed configurations in the same way that a physical packing process does. 
Other deterministic rearrangement procedures use more complicated rearrangement rules 
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which may depend on all neighbours (18 ll ) or on fixing the contact network early in the 



rearrangement (82). The latter strategy can be used to maximize disorder in planar disc 



packings (1831 ) . 

Rearrangement algorithms can also be inspired by models for the motion of gas molecules, 
but these models do not represent the physics of the condensation of a gas. Lubachevsky 



and Stillinger ((84 



85l ) begin with points uniformly distributed within a region. Random 
velocities are assigned to each point, and each point begins to grow a disc at a constant 
rate. When discs collide, the collision is elastic and momentum and energy are conserved. 
After a few thousand iterations of disc expansion, the discs form a near-packing. There is 
a three-dimensional version of this algorithm (jsfj) and a version which can pack ellipsoids 



Rearrangement rules can also be based on the minimization of potential functions. One 
can begin with a uniformly distributed set of sphere centres and then rearrange these centres 
via a conjugate gradient algorithm which minimizes a potential (1881 ) . The potential can 



be one used in a thermodynamic model of a packing which is stable until a yield stress is 
exceeded. 

Rearrangement algorithms can also be written with non-periodic boundaries. They can 

fl 

pack spheres on the surface of a large sphere in R 4 , avoiding any boundary effects (|44l ). 
They can be used as the basis of programs to estimate the most efficient packing of a small 



hard-boundaried region by a fixed number of spheres (1231 ). When combined with elements 



of ballistic algorithms, they can be used to pack cylinders (j89l ). 



3.5 Dynamic Algorithms 

Ballistic and rearrangement algorithms have no basis in the physics of packing formation. 
Physical packings are generally formed by pouring physical objects into a container or into 
a pile. In flow, the objects may be nearly packed or barely interacting. When objects strike 
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confining surfaces or other already-packed objects, they expend energy through erosion, 
through fracturing, through friction, and through deforming themselves and the container. 
Exactly how this energy expenditure happens during physical packing formation is unknown, 
as these processes cannot be observed. 

It is possible to construct models for the formation of packings which incorporate idealized 
versions of the unobservable energy expenditure processes. These models, known as Discrete 



3). 



Element Method (DEM) models, were originally developed for modeling rock mechanics (J9( 
DEM models have been proposed for both piling and container-filling processes for both 



cohesive and cohesionless particles (19 lh . They can be adapted to include processes which 



happen in viscous fluids and in environments where both viscous fluids and air are present 



(I92l). Since they are all based on unverifiable assumptions about interparticle behaviour, it 
cannot be assumed that any of them will capture the physics of packing formation. 

The first DEM packing model was that of Yen and Chaki (1931 ). which simulated the 
formation of a packing in a rigid-walled box. The spheres are initially positioned above the 



box in a configuration produced by a simple sequential inhibition process (I94l ). As time 
evolves, gravity acts upon each particle. When particles come into contact, they experience 
a Herzian force associated with particle deformation and a tangential frictional force. A Van 
Der Waals force can be introduced to provide a cohesive force between the spheres. When 
frictional forces are omitted, the algorithm produces realizations of estimated mean volume 
fraction 0.633 with higher volume fraction variance than expected. When frictional effects 
are added, the mean volume fraction falls to 0.578. 

Since the initial model of Yen and Chaki, many improvements have been made to DEM 



models ((92 



91 



). Rotational interactions have been added ((951), as has the capacity for 



packing polydisperse spheres j96j ) and complex shapes (1971 ). 
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4 Descriptive Statistics 

Any one realization of a packing from a packing process consists of a long list of object loca- 
tions, each of which may be accompanied by size, shape, and orientation information. To fit 
a model to a physical packing process, it is necessary to have statistics which can summa- 
rize what each realization from the same model have in common, and which can distinguish 
between realizations from different models. In physical and engineering applications, it is 
necessary to have statistics which can characterize the physical response of interest and also 
to have predictive statistics which explain how the response depends upon packing structure. 
These statistics may be spatial averages of various kinds, but may also may be measures of 
variability of extreme behaviour. 

In physical examples, the search for descriptive statistics is analogous to the process 
which developed the original thermodynamic state variables used in describing ideal gases. 
Unlike the statistics in classical thermodynamics, descriptive statistics for packings must be 
thought of as random variables having distributional properties beyond their means which 
must be modeled. In the physics literature, it is customary to estimate a statistic from a 
single large sample, and then to equate the estimate to its expected value via an implicit 
appeal to a law of large numbers. This may be reasonable in some contexts, but for any 
one statistic there is no simple way to establish how large the packing must be before the 
between-realization variance of that statistic becomes negligible. 

To be useful for model fitting, descriptive statistics must be able to distinguish between 
realizations from different models in spite of between-realization variability. Statistics for 
spatial processes often lack the power to make such distinctions, and so many statistics may 
have to be tried before a collection can be found that are sufficiently powerful. To find 
these statistics requires having a library of statistics which can describe packings in as many 
different ways as can be conceived of. Statistics are needed to describe what the expert eye 
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can see, and also what no human eye can identify. One of the most significant challenges in 
developing inference for packings is to develop statistics which can be used in these libraries. 



4.1 Random Set Statistics 



Sphere packing processes are examples of random closed sets (J98 



991 ). The random set $ is 



formally defined by the collection of random indicator functions 

I$(x) = 1 if x G $ 
= if x g $ 

for all set 11 . Taking the expected value of this indicator function and its products allows 
moments M p to be defined at every point {xi, . . . , x p } G R d . 



M p (xi, ...,x p ) 



E[I^(x!) . . . h{x p )} 

Pr[xi G $ and x 2 G $ and . . . x p G 



These probabilities are defined on ensembles of many realizations, and do not necessarily 
describe any particular aspects the internal structure of a single realization. The k th moment 
is known as the k— point correlation function or the k— point probability function in statistical 
physics. 

If the random set is stationary, then reduced moments can be defined: 

m p (xi, . . . x p _i) = E[I 9 (0)I 9 (x - 1) . . . I^ixp-i)] 

= Pr[0 G $ and Xi G $ and . . . x p _i G $] 

The reduced first moment mj is referred to as the volume fraction and 7712(7") as the covari- 
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ance. 

If $ is also ergodic, then the reduced moments can be estimated from a single realization: 



mi 

m2.{x) 



I0n A\ 

\A\ 

|0 n An^ x n A a 



\AnA 



x I 



where A is the window of observation, | A\ is the volume of A, and A t = {x e M. d : x + t G A}. 
These estimators are often themselves estimated by sampling at a regular or random arrange- 
ment of points within A. When the random set is non-stationary, mi estimates 



whose interpretability depends upon the spatial variability of Mi(t). 

Both the packing and the closure of its complement can be considered to be random sets, 
so long as the boundary between them is well-behaved and not some sort of fractal. This is 
the case for a sphere packing, and so moments m c k of the realization of the complement of 
the random set can also be defined. Mixed moments can also be defined: 

m 110 (x, y) — Pr[0, x e <E>, y e $ c ] = m 2 (x) — m 3 (x, y). 

While mixed and complementary moment functions can be calculated from the moments of 
the random set, they may be of use if they can illustrate differences between processes more 
accurately than the moments can. 

The most commonly used descriptive statistic for packings is mi. When the random 
closed packing was believed to be a well-defined physical state, mi showed that the earliest 
simulation algorithms were simulating looser packings. Later simulation algorithms were 
judged by their ability to attain estimated volume fractions close to those of a dense random 
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packing (1761) . This statistic cannot be used alone for characterizing packing processes, since 

mj can be used to 



52). 



it reveals nothing about interactions between spheres. Local estimates o 
describe differences in structure within non- stationary sphere packings (llOO 

Random set second moments m 2 (r) are rarely estimated. Instead, point process statistics 
are used to describe sphere interactions. For monodisperse sphere packings, there is a ver y 



close relationship between m 2 (r) and the pair correlation function for the sphere centres (11 Oil ). 

Moment statistics can also be applied to transformed random sets. If the spheres grow 
at a constant rate until they fill space, estimates of and the Euler-Poincare coefficient 
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can be plotted as a function of the degree of expansion (1102 ; 

The complement of the packing is its void structure. This structure can be described 
by the spherical contact distribution function S(r), which is also known as the pore size 



distribution function ((73 



82j). Given an arbitrary point x in the void, S(r) is the probability 
that the nearest point on a sphere to x lies within a distance r of x. It can be estimated 
by finding the distance to the nearest point on a sphere from many locations which have 
been randomly or systematically sampled from within the void. It has been used to compare 
physical packings with simulated packings and to investigate th e app licability of theoretical 



approximations for S(r) which arise from statistical mechanics (Il04l ). 



4.2 Point Process Statistics 

The centres of each sphere in a packing constitute a realization of a point process. The 
basic point process statistics are based on the analogues of the first an d sec ond moment for 



a random measure, and on the Palm distribution of the point process ( 11051 ). 

The intensity A, defined for a stationary process to be the mean number of points per 
unit volume, is seldom estimated since A is a constant multiple of mi. Local estimates of 



intensity can be used to quantify disorder (1361 ) and to investigate the internal structure of 



large physical packings (1291 ) . 
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The K— function K(r) is the second reduced moment function of a stationary ergodic 
point process. The quantity \K(r) is the expected number of points in a disc of radius r 
about a randomly chosen point, which is not counted in the expectation. The if— function is 
seldom used in the study of packings, and the pair correlation function g(r) is used instead. 
It is defined by 

1 dK(r) 

g{r) - 



dbar^ 1 dr ' 

where is the dimension of the unit ball in M d . The pair correlation function may be also 
be referred to as the radial distribution function, although that name is also applied to the 
quantity RDF(r) = \dbrfr d ~ 1 g{r). The if -fu nction has also been generalized to a function 



defined on the number of r— close triples (11061 ) in a realization. Estimation of the if -function 



for stationary ergodic point processes requires careful treatment of window edge effects (11071 ). 

If the point process is not stationary, estimators of #2( r ) for a stationary ergodic process 
are difficult to interpret. Estimates are made by taking each point on the realization and 
creating a sequence of thick shells around it. For each shell, the number of other points in 
the shell is counted. These counts are totaled for each shell type and rescaled based on the 
shell volume. A histogram of the rescaled totals is plotted. If the packed spheres have unit 
radius, a plot of this estimate has a maximum at r = 1, then has a pair of very sharp local 



maxima at at r = 1.73 and r = 2 (1431 ). This type of plot can be used to show the transition 



from a packing to a looser structure in a DEM model as the Van Der Waals force increases 



(1951 ) . For planar packings, the estimated pair correlation function can be used to create a 



pattern that would be seen if X-ray crystallographic methods were applied to a real material 



with the same structure as the realization (1481 ). Those familiar with crystallography may 



find these patterns easier to interpret than estimates of #2 (?")• 

Given an arbitrary location x G M. d , the empty space function H s (r) is the probability 
that the nearest point in a realization to x lies within a sphere of radius r centered at x. For 
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a stationary point process, the nearest-neighbour function D(r) is the probability that the 
nearest neighbour to any point in the realization lies within a distance r of that point. For 
the Poisson process, these two prob abilit ies are identical. For all other point processes, they 



can be compared by the J-function (11081 ) . defined as 

r , x l-D(r) 
J(r) 



1 - HJr) 



Nearest-neighbour functions can also be extended to the 2nd-nearest-neighbour, the 3rd- 
nearest neighbour, and so on (Q). Estimation of all of these quantities require careful 
attention to window edge effects. 

4.3 Statistics Based on Triangulations and Tessellations 

The triangulation of the sphere centres in a packing forms a simple and natural description 
of packing structure. Near neighbours can be clearly defined as spheres whose centres are 
connected by a triangulation edge, and the contact network of the packing is a natural 
subgraph of the triangulation. The dual graph of the triangulation is a simplified description 
of the void structure of the packing. Statistics based on triangulations and tessellations were 
initially applied to packings by physicists who were using phy sical packings as models for 
the molecular structure of liquids and amorphous solids (|8|; l43l ). 

The Delaunay triangulation and the Voronoi tessellation are generally used as bases for 
statistics. The tessellation is constructed by finding all points in R d which are equidistant 
between sphere centres. Euclidean distance is used in the construction, but other distance 



measures can be used to generate different tessellation structures (119l). The triangulation 
is the dual of the tessellation, generated by joining pairs of sphere centres which define a 
tessellation edge. 

The simplest statistics that can be extracted from triangulations and tessellations are lists 
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of characteristics for each cell. For both types of cells, the area, the perimeter, the largest 
and smallest angles, and the longest and shortest sides can be found. The coordination 
number of each sphere is defined to be the number of triangulation edges which belong to 
the contact network. The lists of observations are traditionally summarized by histograms or 
by summary statistics (mean, standard deviation, minimum, and maximum). In the physics 
literature, the histograms are often referred to as plots of the distribution of the statistic. This 
statement is misleading, since the data is a collection of spatially dependent observations. In 
these cases, the histograms summarize aspects of the structure of an individual realization 
rather than estimating a population parameter. For Voronoi c ell a reas, Gamma densities 
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have been found to summarize the shape of the histograms (1831 ; 

Tessellation statistics have been used to investigate the differences between realizations 
of different physical and simulated packings. Differences between packings are generally 
specified first in terms of differences in the estim ated volume fraction of spheres m{. Plots 
can be made of mean contact number versus fru_ 



Voronoi cell face area and volume against m\ (jllll ). Comparisons of histograms of contact 
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29l ). or of the standard deviations of 



numbers from different algorithms are not powerful enough to clearly distinguish realizations 



from different algorithms (I53l ) 



Triangulation-based statistics can also be used to identify possible rattlers. In a physical 
packing, there are often a small number of spheres which are not immobilized. The packing 
surrounding the rattlers consists entirely of jammed spheres forming arch structures, also 
known as a bridgin g structur e s. T hese structures cause blockages when silos filled with bulk 



solids are emptied (1112 
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14). The fraction of rattlers in a packing can be used as a 
measure of packing efficiency (81). Identification of rattlers is complicated by uncertainty 
about whether or not neighbouring spheres are in contact. 

More elaborate statistics can be developed from the tessellation and the triangulation. 



The local density of a packing can be defined as the ratio of the sphere volume to the volume 
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of its Voronoi cell (1291 ). The escape fraction statistic is the empirical cumulative distribution 
function of measurements of the largest sphere that could escape from each tessellation cell 
through the gaps between neighbouring spheres. It has been shown to distinguish between 



physical packings of different volume fraction (J29j). A topological density has been defined 
which is based on a notion of topological distance. This distance is defined by choosing 
a sphere, and then identifying a sequence of shells radiating out from it. The first shell 
contains only those spheres which touch the central sphere, the second contains all sp heres 



in contact with the first shell but not contained within the first shell, and so on (11151 ). The 



density is defined as the leadin g coe fficient of the quadratic fit of the number of spheres in 



each shell to the shell number (11161 ). The smallest values of the density are those of point 
lattices. 



4.4 Statistics Based upon Measures of Order and Disorder 

When physical packings were first proposed as models for the molecular structure of liquids, 
researchers sought to determine whether or not packings possessed some form of local order 
(p). This local order would take the form of small subunits with near-lattice structure, 
combined in some complicated way to produce the general disorder of the packing. 

The presence of an ordered structure in point patterns ma y no t be obvious to the eye. 



118l ). The arrangement 



Materials have been found in nature which are quasi-crystalline ( 1117 ; 
of atoms in these materials would appear disordered to the eye, but x-ray diffraction reveals 
that their structure can be modeled by a projection into R 3 of a higher dimensional regular 
point lattice. 

Disorder can be defined as the absence of a regular point lattice structure for the sphere 
centres. The definition of disorder is not constructive, but instead merely identifies all 
configurations of spheres which are not point lattice structures. While it may be possible 
to define statistics which identify the degree to which a particular configuration is distinct 
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from a regular lattice, there will be many different ways of measuring these distinctions and 
none of them will be clearly better than any other. 

It is necessary to distinguish between topological and geometric disorder. A packing is 
topologically ordered if the sphere centres can be continuously translated s o as to transform 



its Delaunay triangulation into a point lattice without breaking any bonds (Ill9l ). In a plana r 
packing, topological defects of the lattice structure can be easily identified and counted (11201 ). 

A packing is geometrically disordered if any aspect of the arrangement of the spheres 
differs from that of a point lattice. If the spheres in a point lattice are reduced in size by 
a small amount and then randomly perturbed over a small distance, then topological order 
can be maintained in a geometrically disordered packing. Statistics which quantify geometric 
disorder are based on local measurements which describe deviations from expected properties 
of point lattices. 

Statistics can be derived from the locations of contact points on individual spheres, 
expressed in spherical coordinates. The fourth- and sixth-order spherical harmonic functions 
can be evaluated for the contact points on each sphere, and then averaged either over the 
individual spheres or over many spheres. The fourth-order harmonics have non-zero averages 
in the presence of local cubic lattice structure, while the sixth-order harmonics have non- 
zero averages in the presence of local hexagonal close-packed lattice structure. Averages 
over single spheres have been used to study the structure of very large physical packings 



(I29l). Averages over many sphere s we r e ori ginally developed to study the emergence of 



crystallization in simulated liquids (1121 



122f). and improved averages were used to compare 
realizations of simulated sphere packings (1361). Averages over many spheres based on the 
sixth harmonic were used to describe the reduction in disorder observed in simulated sphe re 



packings that had been cycled through the Jodrey-Tory packing simulation algorithm !] 1231 ) . 



Two different many-sphere averages of the si xth harmonic were used to compare realizations 



from three different rearrangement models (1l24l ). Neither average was powerful enough to 
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distinguish between realizations from different models. 

Statistics based on spherical harmonics cannot identify ordered structure that is found 
within small clusters of neighbouring spheres. Statistics which reveal this type of ordered 
structure can be constructed using the side lengths of Delaunay simplices, which are tetra- 



hedra formed by four Delaunay triangles which share common edges (11251 ). In studies of 
the reduction of disorder over time in Jodrey-Tory packings, measures of tetrahedracity and 
quadroctahredracity were found to be more powerful at tracking changes than were were the 
averages of spherical harmonics (1801 ) . 

Other measures of disorder can be found by estimating means and variances of local 



estimates of intensity (1126 



123) 



4.5 Statistics Based on Models for Physical Properties 

Mathematical models for physical phenomena can be applied out of physical context to yield 
new statistics. If a composite of two electrically insulating materials can be simulated by a 
sphere packing, then one or both of its phases can be assumed to be conducting. The bulk 
resistance of the material can be estimated if the internal structure is known. This resistance 
has no physical meaning, but is based on the solution to a set of partial differential equations 
which use the packing structure as part of their boundary conditions. 

Statistics based on models for physical properties are traditionally used as response statis- 
tics in physical applications. When searching for statistics to explain the response, statistics 
based on physical models used out of context can be proposed as potential predictors of the 
response. 

4.5.1 Statistics Based On Prictionless Flow 

In frictionless flow, the flow experiences no internal resistance due to shear. Heat flow by 
means of conductance and the flow of electricity are examples. Models for frictionless flow 
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are relatively simple to construct. 

For electrical flow, the packing and its complement can be considered to be two materials 
with differing electrical conductivity. On the boundary of the packing, two disjoint sets of 
spheres can be considered to be connected to electrodes of infinite conductance. When a 
unit potential is placed between these electrodes, the potential can be calculated at all points 

re c omposite can be 



within the packing and its complement. From this, a bulk resistance for t 
calculated. This resistance defines a mean distance across the packing (11281 ). whose form is 
determined both by the relative resistivity of the two phases and by the disordered structure 
of the packing. 

For a composite modeled by a disordered packing, the bulk resistance is difficult to cal- 
culate. There are no exact methods, and numerical methods require that the packing be 
discretized very crudely. If equal resistance is assigned to all of the edges of the contact 
network, then the potential at all vertices of the circu it can be found easily be means of the 



properties of random walks through the network (I129I ). The bulk resistance is a weighted av- 
erage of potentials at the electrode regions. If the network lies in the plane, the potentials can 
be plotted and used as a diagnostic tool for finding problems with the potential-calculating 
software. A plot of the current along each edge, calculated from potential differences between 
the defining vertices, is more useful than the potential as a descriptor of structure. Bulk 
resistance h as be en used to develop a test for the presence of anisotropy in sphere packing 



realizations (11301 ) 



Models for heat conduction can also be used. The spheres in a packing can be expanded 
in order to generate contact surfaces between neighbouring spheres. A plot of the bulk heat 
conductance of a packing as a function of the degree of spheri cal e xpansion has been used 



to distinguish between packings generated by different models (1l3ll ) 
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4.5.2 Statistics based on shearing flows 



If a fluid flows through a fixed packing, it develops internal frictional losses which depend 
on the structure through which the fluid flows. If the packing itself is made to flow as a 
powder, then its flow is strongly affected by frictional losses arising from colliding spheres. 
These losses greatly complicate the modeling of the flows. 

If the packing is considered to be fixed, the flow of a fluid through its complement can 
be modeled. All fluid flow is assumed to be laminar, since modeling turbulent flow through 
a complex structure is impractical. Major simplifications are required in order to be able 
to calculate bulk flow properties. The complement can be represented by a piping network 
whose structure is determined from the Voronoi tessellation, and pipe resistances can be 



assigned on the basis of local void geometry (1181) . Th e flow through the void structure 



can also be modeled by a lattice-Boltzmann model (jl32j). The velocity profile and pressure 



gradient of the simulated flow can be used as statistics. If the diffusion of contaminants 



through the flow is modeled (11331 ). many different statistics related to contaminant flux and 
concentration can be calculated. 

If the packing itself flows, this flow can take place in a vacuum, in air, in a liquid, or 
in both air and liquid. These flows can be modeled using DEM mo dels which represent 
compac tion processes that do not induce particle fracture 



drums (11361 ; Il37l ; 



1381 ). and fl ow during avalanches (1139 



471). and shear test apparatus (11431 ; 



134 



1351 ). flow in m ixers anc 



14fJ). Flow in mutiaxial (1141 



142 



1441 ) can also be simulated. In some cases, the simulation 
program would impose conditions on the structure of the sample. The best statistics would 
emerge from models for flows which are nearly packed. Flows of fluidized packings would 
be less useful, since the packed structure would almost immediately be destroyed by the 
fluidization process. 

Use of DEM models is complicated by the complexity of modeling nearly packed states in 
motion. In the applications literature, the unverifiable assumptions made in order to obtain 
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simulation results are not often clearly stated. Modeling some aspects of packing formation 
remain primitive. In physical powders, some energy is expended through deformation of the 
grains. Modeling of this deformation can only be underta ken in the plane by means of finite 



146T ). More idealized models 



element methods, and only for a small number of grains (1145k 
can be used in model assessment. A plot of the number of sliding contacts as a function of 
the total number of cont acts during simulated flow can distinguish between realizations of 



different packing models (11471 ). 



4.5.3 Thermodynamic Statistics 

Given a packing of spheres, each sphere can be very slightly shrunken in place. The shrunken 
spheres can be acted upon by a pre-specified potential, and the motion of the particles under 
that potential can be tracked. The motion is confined to a region with periodic boundaries. 
Algorithms for simulating the se motions ha ve been developed as part of thermodynamic 



models for hard-sphere atoms (1148 



149; 



150l ). and also for non-spherical shapes (15 Ol). Once 



these algorithms have been run for some time, statistics representing pressure, temperature, 
and other thermodynamic quantities can be calculated. 



4.6 Graphical Methods 

New statistics can be developed through quantification of characteristic features of patterns 
seen in plots of packings. A simple plot of the discs in a realization can reveal the near-ordered 



arrangements of some planar packings. 



n space, cross-sections showing features of structure 



can be made from tomographic data (l29l ). Two-dimensional projections of tomographic data 



the packing near 



from physical sphere packings in cylinders reveal the ordered structure o: 
the cylinder wall, in contrast to the disordered structure near the centre (1321 ) . 

In the planar case, extra features can be added to images of the packing. The Delaunay 
triangulation can be plotted on top of the discs, and triangles can be coloured to identify 



30 



triples of discs in mutual contact (1861 ). Edges of the triangulation which are associated with 
defects can be coloured to illustrate the overall structure of those defects in large packings 

14- 



4.7 Other Statistics 



There is no evidence that statistics based on presently existing physical and statistical theory 
will be powerful enough to distinguish between the realizations of a given pair of packing 
processes. 

Experimental physics is an important source of new statistics. Size segregation in agitated 
granular matter is a phenomenon that would not be predicted from theory or from simple 
intuition. It occu rs when physical packings of particles of differing sizes or densities are 



shaken (1151 



1521 ) . An example of this is the Brazil nut paradox, in which the large and 



heavy Brazil nuts rise above smaller and lighter nuts when the container of nuts is shaken 



IM). These 



(1421). Simulated sphere packings have been used to investigate the paradox (1153 ; 
simulation experiments could be used as sources of new statistics, particularly in cases where 
there is a size distribution among the spheres. 



4.8 Smoothly Defined and Structurally Defined Statistics 

A statistic is smoothly defined if it is based on measurements taken at many locations in 
the realization, and if the calculation of the statistic treats these measurements as if they 
arose from a random sample. The observation locations can be taken on a fine grid or 
from a realization of a Poisson process, but their definition is completely independent of the 
structure of the packing. To estimate a stationary 2-point correlation function at the vector 
r, a fine grid of sampling locations and its translate by r are imposed upon the realization. 
For each sampling location and its translate, a 1 is assigned if both points are in spheres, 
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else a is assigned. The estimator is the average of the measurements. While this statistic 
is considered to be a simple measure of interaction, its construction combines information 
about voids and neighbouring particles with no consideration of the overall structure of the 
realization. 

Other smoothly defined statistics can be constructed from the Delaunay triangulation or 
the Voronoi tessellation. When a histogram is made of all triangle areas from a triangulation, 
all information about how the triangles fit together is lost. The estimation of the nearest 
neighbour function at a distance r is done by finding the length of the shortest triangulation 
edge from each vertex, and assigning a 1 if it is less than r in length and a otherwise. The 
estimator is again an average, and all other information about structure is ignored. 

Statistics based on physical properties differ from the smoothly averaging statistics in 
that their construction depends on the disordered structure of the entire realization. Bulk 
electri cal r esistance can be thought of as a distance across the network on which current 



flows (11281 ). This distance is based on a weighted average of all the possible paths that could 
be taken across the network, with the weights being chosen to be consistent with the laws 
of electricity and the paths being defined by the disordered structure of the material. This 
would also be true of any statistic based on path properties of a random walk, and also 
on the properties associated with mechanical deformation. The form of these statistics are 
structurally determined by the entire disordered structure of each realization. 

Statistics which are smoothly defined are ideal for use with the Poisson process, in which 
there is no interaction between the observed points. For a point process defined by the centres 
of packed spheres, smoothly defined statistics are often estimated on a length scale which 
is much smaller than that of the individual sphere diameters. Information about individual 
sphere shape becomes confounded with interaction effects. This confounding makes the 
value of the statistic difficult to interpret in terms of features of the realization which can 
be seen. Structurally determined statistics are determined by the disorder of the model. 
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When physics-based statistics are used on packings, all of the information is captured on the 
same length scale as the disorder. There is a vast literature associated with the electrical, 
granular-flow, and fluid-flow models which relates the models to easily-visualized physical 
phenomena. Structurally-dependent statistics may be more useful than smoothly-defined 
statistics at summarizing complex structural details that can be identified by eye. 



4.9 Asymptotics 



Finding the exact distribution of point process statistics is very difficult or impossible, except 
in the case of the Poisson process. If the realizations were very large, then it may be possible 
to find asymptotic between-realization distributions for some statistics. 

The proofs of central limit theorems and laws of large numbers depend on the individual 
observations being weakly dependent. If a point process statistic is a smoothly defined 
average of local measurements, then a central limit theorem can be found as long as the 
pair correlation funct ion for the point process and some similar higher moments decay to 



sufficiently fast (11551 ). This result is independent of the nature of the disorder within the 



realization. If the statistic is not a smoothly defined average, it is not clear that a central 
limit theorem could be proven. 

Laws of large numbers are of great importance in materials science and physics. If a 
composite can be modeled by a packing with the spheres and the void having different re- 
sistivities, then sufficiently large specimens will have the same bulk resistance when a limit 
law for spe cimen resistan ce exists. In the case of ideal gases, equilibrium statistical thermo- 



dynamics (1156 



157 



1581 ) derives limit laws from probabilistic models of atomic behaviour. 
The limiting quantities, experienced as pressure, temperature, volume, and entropy, show 
no signs of variability arising from the dynamic chaos present at the atomic level. Attempts 
have been made to c reate thermodynamic models for the bulk physical properties of granular 



materials (1159 



160), but these attempts have been based on the assumption that any such 
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model would take the same general form as classical thermodynamics. There is no reason to 
believe that limit laws for granular materials would resemble those for ideal gases, since the 
microscale interactions of grains differ from those of atoms and molecules. 

Laws of large numbers are not applicable in many engineering problems. In the case of 
gases, the limits are taken over collections of more than 10 15 atoms. For granular and com- 
posite materials, only thousands or millions of particles will be present, and the interactions 
among those particles will be much stronger than those found in gases. This suggests that 
the proper limiting behaviour is not that of convergence to a number, but rather conver- 
gence to a random variable. In the absence of theoretical arguments for constructing these 
limits, limiting distributions for non-average-based statistics will need to be inferred from 
large samples of model realizations. 

The mathematical intractability of probabilistic limit laws for disordered spatial patterns 
has resulted in the development of alternative approaches to finding constant large-sample 
limits for statistics. If a composite is modeled by a lattice packing of spheres, then methods 



from potential 



conductance (1161 



the ory c an be used to get very good approximations for the bulk electrical 



1621 ). When the packing is disorder ed, t hese methods cannot be used 



1641 ) can be used to find the 



on account of the lack of symmetry. Homogenization (1163 ; 
limiting behaviour of statistics defined by differential equations which have spatially varying 
boundary conditions. Often, these results are presented as bounds between which the limiting 
value of the statistic must fall. While the mathematical arguments behind homogenization 
represent an elegant use of the calculus of variations, the underlying assumptions for these 
arguments generally cannot be related to the microstructure of any real material. In one 
case, a mathematical argument has been constructed to establish that a microstructure exists 
whic h has the coherent potential approximation for bulk electrical conductance as its limit 
( 116.4 



34 



5 Inference 



The central problem of inference for sphere packings is model assessment. No proper statis- 
tical methods yet exist for determining if a fitted model for a packing has anything beyond 
a superficial resemblance to the physical system that it is supposed to represent. This has 
serious implications in physics and engineering, where complex simulation models based 
on unverifiable assumptions are being used. The history of thermodynamics suggests how 
methods for assessing fitted models can be developed. 

No algorithm for the simulation of a packing creates realizations in a manner which 
faithfully represents the physics of packing formation. The DEM models come the closest, 
but even these are based on assumptions about unobservable mechanisms of energy loss 
through friction and deformation. Parameters are few in number relative to the size and 
complexity of the realizations. Parameters often have no physical meaning, especially in the 
case of rearrangement and ballistic models. 

Model parameters can be fi t by methods similar to the method of moments, such as the 



method of minimum contrast ( 11661 ). This method is based on choosing a set of parameter 
values which minimizes the difference between a statistic calculated from the model and 
a statistic calculated from the data. Often, the statistic chosen is the K— function. This 
method ensures that there is some resemblance between the model output and the data, but 
it says nothing about how far that resemblance can be extended. Since the K— fu nction lacks 



the power to distinguish point patterns which can be distinguished by eye ( 11671 ). there is no 
reason to believe that a fit based on it will result in a good fit for other statistics calculated 
from the packing. 

A fitted model must be assessed to determine if the model is useful or if the fit is su- 
perficial. Since every realization from the model or from the data may be different, this 
assessment must be carried out using descriptive statistics which can summarize common 
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aspects of different realizations of the same process and which can identify evidence of impor- 
tant differences between the model and the data. Since there is no guide to which statistics 
are useful, it is necessary to calculate many different statistics for each realization and then 
to compare their distributions by means of multivariate statistical methods. In physical 
applications, the need for many statistics will require the development of accurate methods 
for imaging the three-dimensional structure of static and flowing physical packings. 

Since no model directly and accurately represents the formation of a physical packing, 
it is likely that at least one statistic can be found which gives evidence that realizations 
of a model differ from those of the data. This difference only matters if that statistic has 
an effect on the physical property being modeled. In a typical application, the packing is 
being used as part of a model for some specific physical phenomenon such as the electrical 
conductance of a composite material. If the phenomenon of interest is summarized by a set 
of response statistics, these statistics can be used to fit the parameters of the model. The 
only descriptive statistics that are needed for model assessment are those which affect the 
distribution of the response statistics. In general, there is no way to theoretically determine 
all of the statistics which can affect the response. Further modeling and testing will be 
required to establish if statistics which identify a difference between the data and the model 
also affect the response. 

Some progress has been made towards the development of methods for model assess- 
ment in packings. Various methods have be en de veloped for describing deviations from 



complete spatial randomness in point patterns (1l68l ). but it is not clear if these methods will 



be useful for distinguishing between d iffer ent types of non-Poisson disorder. Residuals for 



point processes have been developed ( 11691 ) which can be used to detect spatial trend and 



interpoint interaction effects. Physicists have used descriptive statistics to quantify differ- 



ences between simulation models for packings (1170 ; 
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1721 ) and between different types 



of physical packings (l29l ). Descriptive statistics have also been used to quantify the effects 
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of changing parameters in DEM models (11731 ). These efforts have been used to suggest that 



some models may be physically wrong, but this condition is too strong. The identification of 
differences between models is only the beginning of the development of objective and useful 
model assessment procedures. 



6 Conclusions 

Sphere packings form a very useful class of models for physical phenomena which cannot be 
modeled formally. They form a class of statistical models consisting entirely of physical and 
computer simulation procedures which need to be fit to data. 

When packings are used as models, the fitted models must be further assessed to deter- 
mine whether or not the resemblance of the model to the physical system is superficial or 
useful. If this is not done, there is a severe risk that the fitted model may be misleading 
as a physical model. The development of assessment procedures requires taking the same 
approach as early physicists; a specific response must be defined, statistics which affect the 
response must be identified, and these statistics must be compared for the model and the 
data. The phenomena being modeled are as well-understood as the properties of gases were 
in the 18 th century, and require a similar experimental approach to identify a useful way to 
summarize everything of importance about them. 
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